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Q ■ A new spectral method is built resorting to (0, 2) Jacobi polynomials. We describe the origin 

O 
(N 



< 



and the properties of these polynomials. This choice of polynomials is motivated by their 
orthogonality properties with the respect to the weight used in spherical geometry. New 
results about Jacobi-Gauss-Lobatto quadratures are proven, leading to a discrete Jacobi 



^ . transform. Numerical tests for Poisson problems in a sphere are presented using the C++ 
» . 

' library LORENE. 
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QQ ' 1. Introduction 
O ' 

. Amongst the various techniques used to discretize partial differential equations, spectral 

o ■ n 

methods, introduced by D. Gottlieb and S. Orszag in the seventies 6|, offer high accuracy 
^ ' at a low computational cost. Their principle is to approximate solutions of PDEs by trun- 



cated Fourier series or high degree polynomials. These methods have an infinite order of 
convergence, because the error between the exact and discrete solutions is only limited by 
the regularity of the exact solution. 

Spectral methods resort to Fourier series in the case of periodic boundary conditions; 
whereas polynomial approximations use orthogonal polynomials (such as Chebyshev or Leg- 
endre polynomials). These methods also employ quadrature formulas to compute integrals 
in variationnal formulations of PDEs. 
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In the present work, we explore a new family of orthogonal polynomials : the (0, 2) 
Jacobi polynomials. After a brief summary about spectral methods in Sec. [21 we detail 
several general properties of the Jacobi polynomials, and present the necessary tools to build 
a spectral method based on the (0, 2) Jacobi polynomials in Sec. O Finally, we show some 
applications to Poisson problems including numerical tests in Sec. HI Concluding remarks 
are given in Sec. [5], and important technical proofs can be found in the appendix (Sec. lAl) . 

2. Spectral methods 

Both numerical analysis and implementation of spectral methods are based on orthogonal 
polynomials, whose major properties are hereby recalled. We shall then emphasize the 
importance of Sturm-Liouville problems in the case of Legendre polynomials, and mention 
results about polynomial approximation and interpolation. Several textbooks cover this 
domain, for instance , or (3|] , and provide comprehensive proofs of the following results. 

2.1. Orthogonal polynomials 

Let A be the interval (—1, 1) and w denote a weight on A, i.e, Wn J^^w{x)x'^dx < oo 
and w > on A. We define 



One can construct a basis of monic orthogonal polynomials by using a Gram-Schmidt 
orthogonalisation process on the basis x", n > : Pq is fixed to 1, then, assuming the P^, 
< m < n — 1 are known, P„ is chosen by 




(1) 



which is a Hilbert space for the scalar product 




(2) 





We recall the following property of orthogonal polynomials: 
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Proposition 1. For all positive integer n, the roots of Pn are real, distinct and strictly 
bounded by —1 and 1. 

In particular, those polynomials do not vanish at 1, we can thus define any family 
of orthogonal polynomials by imposing their value at 1. Besides, the monic orthogonal 
polynomials satisfy the following induction formula. 

Proposition 2. For all integer n > 2, 



X Pn—l\ -Pn—l) \ J-, / \ ||-Pn— 1 



|2 



Pn{x) ={X- ' „— ' ) P„_i(x) - (4) 



p 2 / p 

II / II -'n-2 1 



Finally, a wide class of orthogonal polynomials belongs to singular Sturm-Liouville so- 
lutions, on which spectral discretizations rely. We present the Legendre case {w = 1 and 
= 1), since a general description would be too long. 

Theorem 1. For all integer n > 0, Ln satisfies the following differential equation 

^ ((1 - x^)L'^ix)) + n{n + = 0. (5) 



The A operator, defined by 



A^ = -^{{l-x'y{x)) (6) 



is self-adjoint in L'^{A) by integration by parts. Besides, it is positive and of singular Sturm- 
Liouville type. The equation ([5]) means that Legendre polynomials are eigenvectors of this 
operator, which is the origin of the "spectral" adjective of the numerical methods hereby 
discussed. 

A consequence of ([5]) (through integration by parts) is that, for all positive integers m 
and n : 

L'^{x)L'^{x){l-x^)dx = n{n + l) L^{x)Lr,{x)dx. (7) 



-1 J -I 

This means that the L'^, n > 1 form a basis of orthogonal polynomials in L'^_^2 (A) . 
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2.2. Polynomial approximation error and Sturm- Liouville problem 

This section illustrates why spectral methods are so accurate. The Sturm-Liouville oper- 
ator A enables us to put some upper bounds on the distance between functions with a given 
regularity and a polynomial space, still in the case of Legendre polynomials. This distance 
is computed, by mean of orthogonal projectors. It has been first estimated in [9] and [10|. 

Let be a positive integer, Piv(A) stands for the space of polynomials with degree less 
than A^, and the orthogonal projector from i^^(A) to PAr(A) is denoted by ttjv. For any 
integer m > 0, H"^{A) stands for the Sobolev space of order m on A. 

Theorem 2. For all integer m > 0, there exists a positive constant c, only dependant of m 
such that, for all ip in iJ™(A), 

||V9 - TIn^Wl^A) < cN'^^WipWurr^iA). (8) 

The key of the proof of this theorem is the following one. We denote the n-th 
coefficient of (p in its expansion onto the Legendre basis. Then, 



Since A is self-adjoint in L^(A), 



By iterating k times this result, one can deduce that 

^" = JTT^ , , \,,,M '^\Ln). (11) 

II^"IIl2(A) {n{n + 1)) 

Since (A'^v9|L„)/||L„||^2(yv) is the n-th coefficient of A'^ip in its decomposition over the Leg- 
endre basis, a continuity result on the operator A will allow to conclude. 

2.3. Polynomial interpolation error 

Thus we have seen that any function, given its regularity, can be approximated by poly- 
nomials, with an error decreasing as a power law of the degree of the polynomial. The power 
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depends only on the chosen norm, and the regularity of the function. However, this result 
has few direct numerical applications. Indeed, one has to compute integrals to obtain the 
polynomial approximating the function, which is highly time-expensive. Thus, one shall ap- 
proximate these integrals with quadrature formulas, and replace the orthogonal projectors 
by interpolation operators z^v in the nodes of these quadrature formulas. First, we recall the 
definition of these operators. 

Proposition 3. Given -|- 1 distinct points {xi) in A, and a continuous function ip on A, 
there exists a unique polynomial i^ip in P^v such that 

Wi, iw^ixi) = ^{xi) (12) 

This greatly simplifies the set up of the method, and doesn't reduce the precision ac- 



cording to the following theorem, first demonstrated in 



Theorem 3. For all integer m > 1, there exists a positive constant c, depending only in m, 
such that, for all function ip in if™ (A), we have 

\\^ - iN^Wi^iA) < cN~"'\\(p\\h'^{a). (13) 

The latter result shows spectral methods have an infinite order of accuracy. A C°° 
function is interpolated faster than any power of the discretization parameter. Moreover, it 
can be shown that an exponential convergence of the interpolant is achieved if the function 
is analytic. 



3. Jacobi polynomials 

In this section, several results about Jacobi polynomials are detailed. First, we present 
how Jacobi polynomials arise naturally from Sturm Liouville singular problems. Then, after 
presenting some basic properties, we prove a new result about Jacobi-Gauss-Lobatto quadra- 
tures, enabling us to build a discrete Jacobi transform. Finally, we show some differentiation 
and integration results about (0, 2) Jacobi polynomials. 
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3.1. Introducing J acobi polynomials 

We saw previously that the Liouville operator was crucial to obtain an efficient poly- 
nomial approximation. More generally, a Sturm-Liouville problem consists in looking for 
solutions {u, A) of 



The coefficients p, q and w are three given, real-valued functions such that p is continu- 
ously differentiable, strictly positive in (—1, 1) and continuous at x = ±1; g is continuous, 
non-negative and bounded in (—1, 1); the weight w is continuous and non-negative in (—1, 1). 
One must notice that every Sturm-Liouville problem does not ensure the convergence of the 
associated spectral method. Let us consider the following example ; 



Its eigenvalues are = {j^kY /2, with eigenfuctions (pki^x) = cos((7r/2)A;(x + 1)). Thus, a 
function will be correctly approximated by the cosine series, iff all its even derivatives vanish 
in —1 and in 1. This is due to the non- vanishing of p in the extremities of the interval. This 
Sturm-Liouville problem is regular. Inversely, if p vanishes in —1 and in 1, the problem is 
said to be singular. In the case of singular problems, efficient polynomial approximation is 
ensured. We send back to [7] for a more complete discussion. 

Thus, we are led to consider polynomial solutions to singular Sturm-Liouville problems 
in order to build an efficient spectral method. Let us note (0^) a family of polynomial 
solutions of degree of a singular Sturm-Liouville problem with eigenvalues A^, then one 
can see that q/w is a polynomial of degree zero, p/w has degree 2, and p'/w has degree 1, 
according to the following identity: 




Suitable boundary conditions for u 



1,1) 



(14) 




(15) 



yk (f)k 



XkW 



p 



'k 




<l>'k + 



XkW 



'k- 



(16) 
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Since p vanishes in —1 and in 1, p/w = cq{1 — x){l + x). Let us note ax + b = p' /w, 
then p'/p = {ax + &)/[(! — x)(l + x)] = ai/(l — x) + &i/(l + x). After integration, one has 
p = Ci(l — + x)'^~^^ with ci, a and /5 constants. As a resuh, w = 02(1 — + x)^. 

An integrabihty condition on w imposes a > — 1 and /?>—!. 

We then define the Jacobi polynomials J^'^^ of index [a, (3) as the orthogonal polyno- 
mials for the weight w{x) = (1 — a;)°(l + x)^, normalized by 

" r(i + a)r(n + i)' ^ ' 

where F is the Euler gamma function. 

These polynomials are solutions of the previous singular Sturm-Liouville problem. In- 
deed, ^(1 — a;^)w^ Jn"'^'* j is a polynomial with degree < n, and verifies that for all 
polynomial ip with degree < n — 1 

' ^ ^'1 - x')w^JiC'^A ifw = C 4"''^)--^ ((1 - x^)wip') w = Q. (18) 



l_i w dx \ dx ^ J J_i " w dx 

Therefore there exists A„ such that 



-^({l-x')w^4-'^n=-X^4-'^\ (19) 



w dx \ dx 

Jacobi polynomials are thus the eigenvectors of the previous Sturm-Liouville problem. Iden- 
tifying the highest degree coefficients, one obtains A„ = n{n + a + P + 1). Thus, this 
configuration enables to generalize Theorem [21 

Standard orthogonal polynomials are special cases of Jacobi polynomials. Legendre poly- 
nomials are Jacobi polynomials with index a = (3 = 0. If we denote T„(x) = cos(narccos(x)) 
the n-th Chebyshev polynomial, one can verify that 

= j(-V2,-i/2)r(n + i)r(i/2)/r(n + 1/2). (20) 

Finally, (0, 1) Jacobi polynomials have been studied in Sj. 

3.2. Properties of Jacobi polynomials 

We present here basic results about Jacobi polynomials. It can be proven they are special 



cases of hypergeometric functions (see (l2|]), which provides the following results : 
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Analytical expression 



1=0 



• Induction formula 

2{k + l){k + a + P + l){2k + a + (2:) 

= [{2k + a + (3 + l){a'^ - P"^) + xT{2k + a + (3 + 3) /V{2k + a + j'f''^^ {x) 

-2{k + a){k + p){2k + a + P + 2) j£f (x). (22) 

• Highest degree coefficient 



1 r(2n + a + /3 + l) 
" 2"r(n + l)r(n + a + /3 + l)- ^ ^ 



Norm 



a,m2 2-+^+'T{n + a + l)Tin + P+l) 

n Wlua) (2n + a + /5 + l)r(n + l)r(n + a + /3 + l)" ^ ^ 



• Links between the various Jacobi polynomials families. 

(n + a/2 + (3/2 + 1)(1 - x) 7^"+'''^) = (n + a + 1)7^^^ - (n + . 

(n + «/2 + (3/2 + !)(! + x) 7^^"+') = (n + /5 + 1) J^'') + (n + 

(2n + a + /3) J^"-''^) = (n + a + /3) J^''^ - (n + /3) 4°'?. (25) 

(2n + a + /3) J^''-'^ = (n + a + /3) J^^) + (n + a)jt{^ . 

j("'^)(-a;) = (-l)"J<'^'°)(a;). 

The latter formulae show, for example, that (0, 2) Jacobi polynomials are linked to the 
Legendre polynomials via the following equation : 



{n + 3/2) (1 + x) = (n + 2)L„ + (2n + 3)L„+i + (n + l)L„+2 (26) 
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Figure 1: Plot of J^°'^\ 2 < n < 5. Notice that = (-l)"(n + l)(n + 2)/2 

Finally, as we saw it previously for Legendre polynomials, derivatives of Jacobi polynomials 
are orthogonal for the weight w{x){l — x^), which can also be shown from the identity 

^{jr')=lin + a + P + l)jtV''^'\ (27) 
We show in a very illustrative manner the plot of J^°'^\ 2 < n < 5, in Fig.l. 

3.3. Gauss-Lobatto methods and Discrete Jacobi Transform 

The extrema of orthogonal polynomials are useful in order to construct high precision 
numerical quadrature formulas, i.e, which are exact within a space of high degree polynomials 
: the so-called Gauss-Lobatto formulas. In this section, we recall the definition of these 
quadrature formulas in proposition HI then compute some useful formulas to build numerical 
spectral methods. 

Proposition 4. Let N be a positive integer. Let us denote Xq = — 1 and Xat = 1. There 
exists a unique set of N — 1 points x-i in K, 1 < < iV — 1, and a unique set of N + 1 real 
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Pi, < i < N , such that we have the following identity for all $ in P2Ar-i(A) 



„i N 
J-i 



j=0 

The nodes Xi, 1 < i < N — 1, are the zeros of -^J^n '^^ ■ Moreover, the pj, < j < N are 
positive. 

The quadrature formula 

„1 N 

/ ^{x)w{x)dx c^'^^{xj)pj (29) 

is called Jacobi type Gauss-Lohatto quadrature formula with + 1 points. The expression 
of the coefficients (also called weights) Pj, < j < is given by the following proposition. 

Proposition 5. 

r(iv + i + c.)r(iv + i + /3) 1 



V^n/V-11 2°+/^+^ r(Ar + 1 + a)r(iV + 1 + /3) 1 

^ ' ^' iV(Ar + a + /5 + 1) T{N + l)T{N + I + a + ( (^^))2 " ^^'^ 

n - (.-, I 1 ^ ^"""^^^ r(iv + 1 + a)r(iv + 1 + /?) 1 

^ ^ ^N{N + a + ^+l)V{N + l)V{N + l + a + ^)^jpP)^^^))2- ^ > 

Some commentaries arise naturally at this point. As far as the authors know, this 
expression of the Jacobi Gauss-Lobatto weights has not been given previously in the lit- 
erature. Szego, in gives a very similar expression, but only for Gauss quadrature (i.e, 
the nodes of the quadrature are the zeros of J^)- One can also check the validity of this 
formula for known values of a and (3. For Legendre polynomials (i.e, a = /5 = 0), one 
can check that Vz G [0,A^]pi = 2/[A^(A^ + l)L^(a;j)]. As far as (0,1) Jacobi polynomials 
(denoted M„ hereafter) are concerned, our formulas coincide with those found in 8|, namely 
po = 8/[A^(A^ + 2)M2,(xo)] and Vz G [1,A^]a = ^/[N{N + 2)M%{xi)]. For Chebyshev 

polynomials (« = /? = —1/2), this expression is also valid. Indeed, 
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yiE[i,N-i] Pi 



1 r(A^ + i/2)^ 



1 



mT{N + 1)T{N) (J 



^-1/2,-1/2) I 



(33) 



Then equation ( |20l) allows us to write pj = r(l/2)^/[A^T^(a;i)]. However, in the Chebyshev 
case, Xi = —cos{TTi/N), so that Tjy^Xi) = (—1)^+*. Finally, Wi G [1, — 1] pi = tt/N, and 
Po = pN = vr/2A^. For the numerical applications we have in mind, a = and j3 = 2, so 
that Po = 2A/[N{N + 3)iJp^\xo)f] and G [I, N] pi = 8/[N{N + 3)i/J^'''\x^)f]. 

Now, it is necessary to know how to compute the nodes of the Gauss-Lobatto formulas. 
Looking through the zeros of via a secant method is simple, but zeros tend to accumulate 
near the boundaries of the interval. With no preliminary knowledge of the distribution of 
those knots, it is preferable to use a eigenvalue location method which is enabled by the 
following proposition. 

Proposition 6. The xj , 1 < j < N — 1, are the eigenvalues of the tridiagonal symmetric 
matrix 



with 



7v 



( 



(5i 7i 

7l ^2 

■■. 






















7Ar_2 f^JV-l 



(34) 



(a-/3)(a + /3 + 2) 



(2n + a + /5)(2n + a + /5 + 2) 
/n(n + a + l)(n + /5 + l)(ra + a + /5 + 2) 



1 < n < iV- 1, 



(35) 



1 < n < iV - 2. (36) 



2n + a + /5 + 2\/ (2n + a + /5 + l)(2n + a + /5 + 3) 

The form of the matrix allows for a fast and robust computation of the eigenvalues, e.g, 
by a Givens-Householder algorithm. In the case of a = and /5 = 2, 



— ~/ and 

(n + l)(n + 2)' 



In = 




n(n + l)(n + 3)(n + 4) 
(2n + 3)(2n + 5) 



(37) 
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Finally, these quadrature formulas, enable us to build a discrete Jacobi transform, i.e, 
if one knows the values of a function on the Gauss-Lobatto nodes, on can compute the 
coefficients of its polynomial interpolant in the Gauss-Lobatto nodes. 

Proposition 7. Let us note ijyf = X^fcLo fk^k the interpolant of f on the Gauss-Lobatto 
nodes, then 

Vm < iV - 1, 

2m + a + p + l r(m + l)r(A^ + 1 + a)T{N + 1 + P)T{m + a + p + 1) 



N{N + a + p + l) T{N + l)r(m + 1 + a)T{m + 1 + P)T{N + a + p + 1) 

{N + a + f3 + l) y Jn{xo) JN{xi) Jn{xn) ^ 

Various cases of these formulae can be checked. For Legendre polynomials, they reduce 

to 




iV(iV + l)Z. L%{x.) ' (Ar + i)Z.L;v(a:.)- ^^'^ 

In our case (a = and [3 = 2), 



~ ^ 2m + 1 ^ f{xi)J„,{xi) ~ ^ 1 f{xi) 

iV(iV + 3)^ J^(x.) ' (Ar + 3)^Mx,)- ^^'^ 

This discrete Jacobi transform is an essential ingredient in order to build a Tau method 
a. Indeed, a,. W.on / . .ep.e.,.ed by .he coeffice... /. fro. .he ,a«e. p.opo.«on. 
Differential operators on / can be translated as linear operators on the {fk)-, boundary 
conditions as linear constraints, etc. 

Since we have detailed the properties of Jacobi-Gauss-Lobatto quadrature formulas, we 
may investigate the po^ 



ynomial interpolation error. Theorem [3] has been proven for Legendre 
in Chebyshev in [l3], (0, 1) Jacobi in js], and Gegenbauer polynomials (a = /3) in Q]. 



It still remains to be computed in the (0, 2) case. As an illustration of this section, we 

show the plot of Ji^'^"* and its Gauss-Lobatto nodes in Fig.2, we also include the location of 

Chebyshev Gauss-Lobatto nodes of the same order for comparison. 
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— ito;2) 

^12 

* Jacobi collocation points 
o Chebyshev collocation points 




Figure 2: Plot of J^2'^\ location of Jacobi and Chebyshev Gauss-Lobatto nodes. 



3.4- Inverse inequalities 

The following result is enclosed here for completeness. It may not be of interest for 
our numerical purposes, but may be useful in order to devise some a priori inequalities in 
non-linear problems. See for more details. 

Theorem 4. There exists a real c such that the following inequality is satisfied for all positive 
integer N and all polynomial (fN in P7v(A). 



(42) 



This inequality is optimal. Indeed, there exists a constant c' such that 



\J'n\hI{A) > c'N \\Jn\\lI{A)- 



(43) 
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3.5. Physical motivation for (0, 2) Jacobi polynomials 

Let us have a closer look at the orthogonality property of (0, 2) Jacobi polynomials, 
which we shall denote Jk hereafter. 



Wn,k, J Jkix)Jn{x){l + xYdx = 5k,n\\Jk\\l- (44) 
By mean of the change of variable r = (1 + x)/2, we obtain 

1 

Wn,k, / Jk{2r-l)Jni2r-iydr = -5k,n\\Jk\\l- (45) 
Jo ^ 

Thus appears the crucial weight in spherical coordinates. The utility of such a family of 
polynomial appears when computing global quantities. Let us take the example of a magnetic 
field in a spherically symetric geometry. Its total energy is proportionnal to J B'^r'^dr. 
Denoting 6„ its coefficients in a decomposition on a Jacobi polynomial basis, its total energy 
will then be proportionnal to ^ &^||^n||^- Using Chebyshev polynomials would have required 
to employ desialasing techniques to compute this integral (See j3]). 

Besides, when dealing with non-linear terms, one can use, for example, truncation tech- 
niques, and in our case, truncating Jacobi coefficients doesn't increase the energy of the 
field. For instance, if the previous magnetic field obeys an induction equation with dissi- 
pation, exact solutions have decreasing energy, and we ensure the desaliasing techniques do 
not increase the energy of the numerical solution. 

3.6. Operator matrices for (0, 2) Jacobi polynomials 

Various approaches can be examined when building a spectral method. A C++ library 
called LORENE [l^ has been built and mostly uses the Lanczos-Tau method exposed in [t]. 
Therefore, we need operator matrices such as derivation, integration, multiplication and 
division by X + 1, whose expressions in the (0, 2) case are given in the following proposition 

Proposition 8. For all positive integer n, the following identities are true : 

• Derivation : 
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n-1 



j=o 



^ ' (n + l)(n + 2)_ 



J. 



(46) 



Integration : 



_n+3 1 n 

" " (n + 2)(2n + 3) ~ (n + l)(n + 2)^" ~ (n + l)(2n + 3) ^ ^ 



Division by 1 + x 



l + x 



n-1 



-i-,2j + 3 



(^ + l)(n + 2) (j + l)(j + 2) - 
(j + l)(j + 2) (n + l)(n + 2)_ 



J. 



(4J 



Multiplication by 1 + x 



(I I _ (n + l)(n + 3) I n^ + 3n + 3 n(n + 2) 

^ (n + 2)(2n + 3) (n+l)(n + 2)^"+ (n + l)(2n + 3) ^ ^ 

4. Numerical tests and implementation 

Numerical tests have been performed using the C++ library lorene which pro- 
vides useful tools for numerical relativity, especially elliptic solvers. Poisson equations are of 
outmost importance when looking for solutions of Einstein equations. Indeed, various for- 
mulations of these equations lead to elliptic equations which can be solved through iterative 
Poisson resolutions. 

LORENE performs 3D multi-domain resolutions with spectral accuracy thanks to Cheby- 
shev polynomials in the radial direction, and spherical harmonics or trigonometric poly- 
nomials in the 6 and (f directions. Moreover, it uses several spherical domains, a nucleus 



centered around r = 0, some shells, and an external compactifiec 
a variable). Some applications can be found in Q, Q, and Q 



domain (using m = 1/r as 
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Figure 3: Decay of the errors (maximum difference with theoretical solution on collocation points) in each 
domain. Settings are Nq = 17, — 16. 

We implemented the use of the (0, 2) Jacobi polynomials in the nucleus, where their 
orthogonality properties are the most interesting, and developped a Poisson solver thanks 
to these polynomials. 

The Poisson equation reads 

A/ = S, (50) 

where A = 5^ + ^dr + l/r'^{dg + f^ffi^e + d'^) is the Laplace operator. Our procedure is the 
following one : we perform a decomposition of functions on spherical harmonics, through 

f{r,e,^) = J2fiUr)Yr{0,^), (51) 

l,m 

where Yj^ denote the spherical harmonics. Since AY/^ = + 1)1^™", we reduce equation 
flHJj) to a set of ID equations 



dr'^ 



+ - 



r dr 



-fi 



Im 



Sim- 



(52) 



Then, assuming fim{f) = J2k=o fkimJki^), one can assemble the matrix of the operator 

appearing in (l52l) . and adding some boundary conditions, invert it to find the coefficients 
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Figure 4: Decay of the errors (maximum difference with theoretical solution on coUocation points) in the 
nucleus for Jacobi and Chebyshev methods. Settings are Ng — 17, ~ 16. 

fkim- More precisely, we decompose in three domains, a sphere of radius 1, a shell between 
radii 1 and 2, and an external zone r > 2. We map those three domains to [—1, 1] via the 
following equations : 

for r < 1, r = (l + x)/2, 

for l<r<2, r = (3 + x)/2, (53) 

for 2 < r, r = A/{l — x). 

Then, for each /, we assemble a three- block matrix containing the matrices of the following 
operator in the (0,2) Jacobi basis for the nucleus thanks to ( H6l) and ( l48l) . and in the 
Chebyshev basis for the two other domains 



'^^^[f-/'(-i)]-#±^[/-/(-i; 



:i+^)/'(-i)]- 



(54) 



dx^ ' l + x''dx '' ^ " (1 + x) 

In the nucleus, it corresponds to the operator in ( l52i) . under the assumption that for / = 0, 

/'(-I) = 0, for / = 1, /(-I) = 0, and VZ > 2, /(-I) = /'(-I) = 0. In order to enforce 

those conditions when solving the problem, we replace the last lines of the first block-matrix 

17 



by the values of J„(— 1) or J!^{—1), and those of the source by 0, depending of the value of 
/. Finally, we include matching conditions between the domains in the 3-block matrix. 

Several numerical tests have been passed. First, we used the following smooth source in 
Eq.dnO]) 

5(r, e, if) = 4(r2 - 2 + 3z^)e-''-'\ (55) 
which admits the following solution 

f{r,e,ip) = e-^'-^\ (56) 

We compute the error between the analytical and numerical solution by looking at the 
maximum difference on the collocation points in each domain, which gives an estimation of 
the corresponding norm of the difference. Results are displayed in Fig. 3 and show as 
expected an exponential decay of the error as a function of the number of radial collocations 
points. 

One can compare the efficiency of Chebyshev and Jacobi polynomials for the solution 
(|55l)- (!56l) . When lorene uses Chebyshev base in the nucleus, the mapping is r = Rx, 
where R is the radius of the nucleus, and we use the parity of fim with respect to I to use 
odd or even Chebyshev bases. Notice this isn't possible in the Jacobi case because Jacobi 
polynomials with a ^ (3 do not have a definite parity. Results are displayed in Fig. 4. One 
can see no loss or gain of precision is acquired with this method. In fact, Jacobi polynomials 
are expected to be more useful in dynamical evolutions to ensure stability (see Sec. 13.51) . 

Finally, we examined a non-smooth solution, namely 

S' = 35v^/4 for r < R, S = otherwise 

" (57) 
/ = r5/2 - 7/2^/72 for r<R, f = -hR^'^/2r otherwise 

where R is the outer radius of the shell. 

Results are displayed in Fig. 5, and one can see as expected a geometric decay of 
the error as a function of the number of radial coefficients. We found that the rate of 
convergence in the Chebyshev case is approximately 2.76, and 4.62 in the Jacobi case. 
Indeed, / G Hl^^^_^,^_,^,, and / G 
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Figure 5: Decay of the errors (maximum difference with theoretical solution on collocation points) in the 
nucleus for Jacobi and Chebyshev methods. Settings are Nq = 17, N^, = 16. 

5. Concluding remcirks 

We have described a new numerical spectral method and showed some applications to 
Poisson problems in a sphere. In order to build this method, we have first recalled ba- 
sic results about orthogonal polynomials and generalized quadrature formulas and discrete 
polynomial tranforms. Besides, we computed various differential and algebraic operators. 
Through numerical tests of Poisson equations, we have shown that the method is conver- 
gent and accurate. In particular, for smooth solutions, we find as expected a error decaying 
exponentially with the number of radial coefficients. This method has also been integrated 
to a 3D multi-domain spectral solver. However, the interpolation polynomial error on the 
Gauss-Lobatto nodes remains to be computed. Studying the nodes of the Gauss Lobatto 
quadrature formulas shows evidence that their distribution is very similar to the Legendre 
case, or the (0, 1) case, and therefore, that a similar theorem could be proved. Other appli- 
cations for (0, 2) Jacobi polynomials could be foreseen. For instance, hyperbolic equations 

in a sphere, with applications to fiuid dynamics in a star. 
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A. Appendix 

In this section we provide proofs for the major new results listed in this paper. 
A.l. Proposition [3] 



Proof. To simplify the expressions, J„ will be denoted Jn"'^"* in t 



lis proof, except when the 



preliminary lemma. We 



index is explicitly written. The proof resorts to a already known 
shall note J^^^ = (A„x + ^n)Jn ~ t^nJ'n-i the induction formula of the J^, since they are 
orthogonal polynomials with respect to the weight w{x){l — x^). 
Let us first compute p,, for 1 < i < — 1. 

Lemma 1. 

\ II T l|2 

An-l|K„_i||^2 (^) 

VzG[i,iV-i] P-^ = TT^^rfiiTYf^y (5^) 

Lemma proof . Let us fix i in [1,A^ — 1], et define = ^^^^(1 — x^) in P7v(A). The 

exactness of the quadrature formula ( 128|) then gives 

(1 - x^i)pMxi) = f :^kM(i _ x^)w{x)dx. (59) 

The computation of the integral may be done through the study of the following quantity 

S„(C.,) = ^^±lKl^M^^±lM^. (60) 

C-v 

The induction formula of the allows us to write 

Sn{C,V) = KJniC)J'M + ^nSn-l{C,V)- (61) 
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We here recall that z/„ = A„|| J^||? fA-i/'^n-ill'^n-iIlL fAi' according to the induction 
formula of orthogonal polynomials. Then 



\i\\J'n\\\2 (^) 11-^411^2 (A) -^n-lll-'^-lllia (^) 

Since So{(,ri) = 0, one can see that 



(62) 



Tjj(l — a:^ ) 



This is already known Christoffel-Darboux formula (see 12]) for orthogonal polynomials. 
Let us take that equality inn = A^ — 1, ( = x, r] = Xi, multiply by (1 — x^), then integrate 
with respect to w. Left member is given by 

J'^{x)J'^_'l^{Xi) — J'j^{Xi)J'^_^{x) 2\ / \ T T' f \ ^Af(^) /I 2\ / \ T 

(1 — X )w[x)dx = JM_Axi) / ( 1 — X )w[x)dx. 

1 J —1 ^ 

(64) 

In the right member, by orthogonality, every term in the sum vanishes, except for /c = 1, to 
give 



\n-i II .;nV ' / J[{x){l - X )w{x)dx = \n-i\\J'n-i\\l^ fA) (65) 

IkiIIl2 (A) ^-1 ^{i-x^) 

since J[ is a constant. Then, one deduces that 



1 7/ 



AAr-l||-^^_ll 



^^(1 - X^)wix)dx = ^gQ) 

X Xi ]\j—i\Xi) 

□ 



Let us now compute the various quantities in the lemma. 



Computation of ||t/^_i|||2 



j„(l„^2)(A) 



After integration by parts, and using the differential equation (fT^ . one has 



1 .1 

{JN_if{l-x'^)w{x)dx={N-l){N + a + f3) Jli_^{x)w{x)dx (67) 

21 -^-^ 



I.e, 

Computation of \n-i 

The identification of highest degree coefficients leads to 

^ NkN ^ i2N + a + P)i2N-l + a + P) 

{N-l)kN-i 2{N-l){N + a + f3) ' ^ ^ 

Computation of (1 — x1)J'l^{xi) 

Developing the differential equation ( fT9l) . and evaluating in Xi results in 

(1 - xl)J';,{x{) = -N{N + a + l3 + l)JN{xi). (70) 

Computation of J'i<!^i{xi) 

Induction formula of evaluated in Xi gives 

J'N+l{Xi) = -l^N^N-li^i)- (71) 

Besides, deriving the induction formula ( |22l) of the J„, and evaluating in Xi gives 



, , {2N + l + a + p){2N + 2 + a + p) 
^+1^"^^^" 2iN+l){N + l + a + P) "^^^^^ 



^N + a){N + (3){2N + 2 + a + P) 

JN^iixi)- (72) 



{N + l){N + l + a + (3){2N + a + (3) 

But 2 

^A-II'^^IIL(,_.2)(A) _ (iV + a)(iV + /j)(2iV + 2 + a + /?) 
"""^ A^_i||J]v_i||i WA) iV(iV + « + /?)(2iV + « + /?) ■ ^ ' 

After substitution, one has 
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{N + a)iN + f3){2N + 2 + a + P) 
2N + a + p 



{N + l){N + a + f3 + l) N{N + a + p) 

_ {2N + 2 + a + P){2N + l + a + p) 
2{N + l){N + l + a + f3) 



JN{x^). (74) 



So 



Collecting all above formulas, one has 

_ 2°+^+^ r(Ar + 1 + a)T{N + l + (3) 1 

N{N + a + (3 + l)V{N + l)V{N + a + (3 + l){JN{xi)y ^ ^ 

Let us now compute pN- We apply the quadrature formula (!28i) to $(x) = J^(a;)(l + a;). 
We then obtain 

2p^^(l) = j;,(x)(l + x)u;(x)cix= ^^"^^^^ J^^^l^''^^\x){l + x)w{x)dx. 

(77) 

But (2Ar + 2 + « + /?) J^"'^+'^ = (Ar + 2 + « + /?) 7^,"+''^+'^ - (Ar + /?+ l)4"+''^+'\ according 
to ( l25l) . Multiplying this by (1 + x)w{x) and integrating leads to 

{N+a+(3+2) j' fj^^^'^^^\x){l+x)w{x)dx = {N+P+1) J^^^l'^^^\l+x)w{x)dx. (78) 

Hence 

(79) 

Moreover, 
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Finally, 



Then, 



^ r(a + 2)r(a + 1)T{N)T{N + 1 + (3) 

P""' T{N + l + a)T{N + 2 + a + (3) ' ^ ' 



Since Jiv(l) = r(iiajr("!v+i) ' 

= r 2°+/^+! r(iv + i + a)r(iv + i + /?) i 

^ ^ V(iv + a + /? + i)r(iv + i)r(iv + i + « + /?) (j(;«./^)(3;„))2- ^ ^ 



The expression of po is obtained by exchanging a and /5. 



□ 



74.^. PropositionlE 

Proof. First, let us recall that £ (^jf^^^ = |(n + a + /5 + then we will look 

for the zeros of J^^i'^'^^^ ■ But the J^~^^'^^^^ satisfy the following induction formula 

2(n + l)(n + « + /? + 3)(2n + « + /5 + 2)4°"^''^^'^ = 

[{2n + a + l3 + 3){a-(3){a + (3 + 2)+ xT{2n + a + (3 + 5)/T{2n + a + (3 + 2)] 4"+^'^+^) 

-2{n + a + l)in + p + l)i2n + a + p + A)jtf'^^^^ (84) 

Let us then define 



J. Jr^''^" I i2n + a + p + 3)T{n + l)r(n + a + /? + 3) 

" l|ji"'"''^"''^llL^ (A) V 2"+/^+3r(^ + a + 2)r(n + /3 + 3) 



(1 — a;^)ii? 



Then, the reccurence formula of the J* can be written as such 



J. {a-(3){a + (3 + 2) 

{2n + a + (3){2n + a + f3 + 2) "'^ 



^5) 



2 + a + l)(n + /3 + l)(n + Q + + 2) 

^2n + a + /? + 2\/ {2n + a + p + l){2n + a + P + 3) 



2n + a + /3y (2n + a + /?- l)(2n + « + /? + 1) ""^ ^ ^ 
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or similarly 



which can be written in a matrix form 



X 



Jl{x) 



\ 



( 



5i 7i 

1\ h '■■ 

■•• ■•• ■•• 

■•• 7Af-2 

7Ar_2 5Ar_i j 









Jl{x) 



( 



+ 77V- 








^7) 



Otherly said, the Xj, 1 < j < — 1, zeros of are the eigenvalues of the fl34|) matrix. 



□ 



A . 3. Proposition 

Proof. One can write, for a fixed m < N, 



N 



y^^PiJm{Xi)iNf{Xi 



i=0 



N r N 

k=0 L i=0 



^9) 



For m < N — 1, since k + m < 2N — 1, only /m||<^m||/,2 (^-j is left in the right-hand side. 

Then fm = n, ,X J2iLo Pi'^mixi) f (xi) . Eg. (12^ and the expressions of pi from proposition 

|5]lead to the desired result. 

For m = N, the exactness of the quadrature formula allows us to suppress the — 1 first 
terms of the sum, but the last one is not equal to the norm of Jn- Precisely, the right side 
member is equal to Yld=Q Pi{J N{xi)Y ■ Using the expressions of the weights, one finds 



2^+/3+i Y{N +l + a)V{N +1 + [3) 
■'^N{N + a + (3 + l) r{N + l)r(A^ + ! + « + /?) 

One deduces 



+ 1) + (AT -1) + (« + !)]. (90) 



f 2°+^+^ r(Ar + 1 + a)T{N + 1 + /3) ^ 

N T{N+l)T{N + l + a + [3) ~ /^^^^^^^^^ 



which is the desired result. 



(91) 
□ 



25 



A. 4- Theorem^ 

Proof. First let us compute ||^/j||^2 thanks to the Gauss-Lobatto quadrature with n + 1 
points {J'j^ has indeed degree 2n — 2). 



{n + a + P + iy ( /r(n + /? + l)V fT{n + a + l) 



r(n)r(/? + 2) 



ia+/3-l 



n{n + a + 



1 



T{n)T{a + 2) 
1 \ r(n + 1 + a)r(n + 1 + /3) 



a + 1 13 + 1 J T{n)T{n + l + a + (3) 



Then, for ip^ = '^^=0^"''^^, one has, thanks to Cauchy-Schwarz inequahty 



N 



N 



l<^iv| 



< 



< 



T.\v'f\\Jn\\h 



(A) 



1/2 / TV I 7 |2 

ll'^"llL2,(A)^ 



E 

,n=0 



I I Hi (A) 



1/2 



But, according to the previous computation. 



(92) 



(93) 



Hi,{A) 



^ ^ +^J-^]n{n + a + p+l){2n + a + P + l). 



One can deduce that 



\fN\Hl(A) < IIV^IIl^CA) 



N' 



^ ' n=0 



(94) 



nn + a + /5 + 12n + a + /3 + l 



ll/2 



AT 



< 



1 / 1 



+ 



4 Itt + l 



SUPo<n<Ar 



n n + a + [3 + I2n + a + (3 + I 



N N 



N 



1/2 



(95) 



Therefore, we have 

|<^7vUi(A) < cA^Iv^tvIIl^CA). (96) 

Now, we turn to the optimahty of this inequahty. We will compute ||^Ar||L2_(A) thanks to 
the quadrature formula with A^ + 1 points. 



Il2{A) 



> pM{-^)f + PN{m)f 
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(97) 



because the weights pi are non- negative. Using ( |271) and ( fT?l) . one can find that 



WJ'nWIua) > T^^'^N{N - lf{N + a + l3 + 1){N + « + /? + 2)^ 

T{N + l + a)T{N + l + P) [ 1 , 1 , 



r(Ar + l)r(Ar+l + a + /5) V(/? + l)(/5 + 2)2 (a + l)(a + 2) 
Thus, 

Therefore, 

i^tl^ > c'iV^ (100) 

II "^TN 



'iVllL2(A) 

^ . 5. Proposition 
Proof. • Derivation : 
We recall that 

{2n + 4) J^°'=^) = (n + - (n + 3)4'_:? 

(2n + 3)J^°'^) = (n + 3)Jr)+njf? 
From these equalities, we obtain 

(n + 4)jP) = 2ELo(^ + 2)4°''^ 

{n + 3)(n + 2)(n + 1)7^=^^ = (-1)" ELo(-l)'(2^ + 3)(A; + 1)(A; + 2) jf ''^ 
which allows us to write 

\k 



(n + 4) J(^'3) ^ 2 J](2, + 3)(, + + 2)(-iy jj^ ^^g^j^ ^ [ J^- (103) 



□ 



(101) 



(102) 



The term in between brackets can be computed using a simple element decomposition 
and leads to 

;f (-1)" (-1)" I (-ly (104) 

f-(t + 3)(t+l) (n + 2)(n + 3) ^ (i + l)0 +2)- ^ ' 
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One can deduce that 



+ 4)4''" ^ E(2, + 3) {l + ,-ir"ii±i|±|} jr>, (105) 



Which gives the desired result, thanks to ■£:Jn''^^ = ^{n + 3)J^_f^i . 

Integration : 
Let us recall that 

(2n + 1) J^'''^ = {n + 1) J^') - (n + 1) jf'^ 
2(n + l)J<°'^) = (n + 2)J<°'^)+n4°:¥ 

Then, 

^(^ + 1)4°-? = ^Jt''' = ^^(4''' - 4°:?) 

2 " ^ (is 2A; + 1 (ix " " ^ 

k+1 d f{k+ 2)jr + kjf^^ (k + 1) j£? +ik- 



(106) 



2k + ldx \ 2{k + 1) 2k 

Or 



(107) 



.(0,2) _ f k + 2 (0,2) _ 1 t(0,2) _ k - 1 (o,2)\ . ^ 

'^-i " rfx V('^ + 1)(2A;+1) ' kik + 1) ki2k+iy'-^J- ^ ^ 

By taking n = A; — 1, we obtain the primitive given by the propostion [HI besides, since 
Jn''^\l) = 1, this primitive vanishes at 1. 

Division by 1 + x : 

We will use the following formula : 

(n + 2)(1 + x)4°'3) = in + 3)4°-'^ + (n + 1) ^^J, (109) 

from which 
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But we saw previously that 

{n + 3)(n + 2){n + = (-1)" j^i-lfi^k + 3)(A; + 1)(^ + 2)jf^'\ (111) 



fe=0 

So 

r(0,2) 



- + 3)(n + 2) ^(-l)"-^(2j + 3)(j + 2)(j + 1) 



[n ■ 



1 + x 



V ^ i + '^'^+'^^~-^-' rii2) 

^ (A; + l)2(A; + 2)(A; + 3)2 ( ^' ^ 1 + x ' ^ ^ 



The term in between brackets can be computed using the following simple element 
decomposition: 

4 2 112 1,, 

+ 71^-7X^ + 1^-1^-71^^- (113) 



(A; + l)2(A; + 2)(A; + 3)2 A; + l (A; + l)2 A; + 2 A; + 3 (A; + 3)2 
The sum presents lots of cancellations and is equal to 

4V ^ = ^ ^ . (114) 

{^.{k + lY{k + 2){k + ZY (j + l)2(j + 2)2 (n + 2)2(n + 3)2 ^ > 

re— J 

Finally, 

JSS ^ f (n + 3)(n + 2) _ (j + 2)(i + l) ) (0,2) J^:S{-^) 

1 + x ^ 4 \ (j + l)(j + 2) (n + 2)(n + 3)J ^' ^ 1 + x ' 

(115) 

Notice that this proof can also be achieved using a Christoffel-Darboux formula. 

• Multiphcation by 1 + x : 

We recall that 

(n + 3/2)(l + x) J^') = (n + 2) J^^''^ + (n + 1)^°;? 
2(n + l)J<°'^) = (n + 2)J<°'^)+njS 

Then, 

{n + 3/2)(l + .) J^o.^) = (n + 2) f + '^^^f' ^ + "^-'^ 



(116) 



(n + 3)4|+^+i)^--Y (iir) 
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Or 



X 



,2) 



(n + l)(n + 3) (0,2) + 3n + 3 

{n + 2)(2n + 3) ^ (n + l)(n + 2; 



n(n + 2) 



> + l)(2n + 3) 



(0,2) 



n-1 



:ii8) 
□ 
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